Quantum dynamics in ultra-cold atomic physics 
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We review recent developments in the theory of quantum dynamics in ultra-cold atomic physics, 
including exact techniques, but focusing on methods based on phase-space mappings that are appli- 
cable when the complexity becomes exponentially large. These phase-space representations include 
the truncated Wigner, positive-P and general Gaussian operator representations which can treat 
both bosons and fermions. These phase-space methods include both traditional approaches using a 
phase-space of classical dimension, and more recent methods that use a non-classical phase-space of 
increased dimensionality. Examples used include quantum EPR entanglement of a four-mode BEG, 
time-reversal tests of dephasing in single-mode traps, BEG quantum collisions with up to lO'^ modes 
and 10^ interacting particles, quantum interferometry in a multi-mode trap with nonlinear absorp- 
tion, and the theory of quantum entropy in phase-space. We also treat the approach of variational 
optimization of the sampling error, giving an elementary example of a nonlinear oscillator. 



I. INTRODUCTION 

Quantum dynamics is one of the most fundamental 
problems in modern physics. This is because time- 
evolution is the basis for any theoretical prediction. Yet 
many-body complexity makes this an extremely challeng- 
ing task in quantum systems. New theoretical meth- 
ods are needed, and quantitative experiments with well- 
understood interactions are vitally important in order to 
test predictions. In this article, we review some recent 
developments relevant to ultra-cold atomic physics. 

Ultra-cold atoms provide an exceptionally simple and 
well-understood physical environment, allowing quanti- 
tative tests of dynamical theoretical predictions [l|, 3. 
Recent experiments explore temperatures below InK 3| , 
capable of demonstrating dynamical behavior in many- 
body systems in new regimes. The important new feature 
of these systems is that they allow isolated, macroscopic 
quantum systems to evolve almost unitarily, with very 
little coupling to external reservoirs. It is this feature of 
these experiments which is highly unique, and not found 
in most previous condensed matter experiments 

Features of recent experiments include [5]: 

• Bose-Einstein condensates: atom 'photons' 

• Atom lasers, atomic diffraction, interferometers.. 

• Quantum superfluid fermions: atom 'electrons' 

• Universality: Strongly interacting fermions 

• Superchemistry: Ultracold molecule formation 

• Squeezed BEC: Spin-squeezing with spinor atoms 

An important development is the growing ability of ex- 
perimentalists to measure atomic correlations [6] and per- 
form atom counting experiments with noise levels below 
the standard quantum limit of Poissonian fluctuations. 
A typical schematic picture is shown in Fig ([T]), which 
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Figure 1. Schematic diagram of atom counting experiments 
using metastable Helium and multi-channel plate counters. 



shows a magnetically trapped ultracold atomic cloud. Af- 
ter a dynamical quantum collision of two Bose conden- 
sates, the trap is turned off and atoms are counted by 
the multi-channel plate (MCP) [7] below the trap. 

As well as these atomic correlation experiments, other 
experiments of interest include quantum collision [8l| and 
quantum interferometry experiments. In all these cases, 
there is an external Hamiltonian which can be changed 
non-adiabatically, leading to quantum dynamical evolu- 
tion in the many-body system. This is obtained by ex- 
ternally control of laser or magnetic fields. 

Unlike traditional condensed matter environments, 
these experiments are carried out in a high vacuum, us- 
ing optical or magnetic trapping potentials. It is this 
feature that allows these systems to evolve with almost 
no contact with a heat reservoir. In summary, for the 
first time in physics, we have large many-body quantum 
systems capable of unitary evolution with a wide variety 
of controlled interactions. This creates an unrivaled op- 
portunity for testing calculations of quantum dynamics. 
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II. GENERAL HAMILTONIAN 



In the uniform case, the hopping matrix LUnn' is 



The Hamihonian of the relevant uhra-cold atomic sys- 
tems usuaUy are rather simple, being comprised of well- 
defined single-particle and interaction terms. Thus, 



H = Ho + Hi 



(2.1) 



where Hq and Hj are the non-interacting and interacting 
parts of the Hamiltonian respectively, so that Hq is a 
general linear Hamiltonian, given by: 



2m. 



■^s' (r) dPv . 



(2.2) 

where ^'^ (r) is a quantum field operator r with internal 
spin or atomic species index s — 1, . . . S*, where : . . . : 
indicates normal ordering, and we use the Einstein sum- 
mation convention for repeated indices Q. In addition, 
nis is the mass of species s, Ks'r is a local potential, and 
Hi describes particle-particle interactions with interac- 
tion potential Uss'rr' ■ 

ss' 

(2.3) 

It is convenient to introduce local mode operators to 
treat such quantum field equations [lOj. We describe this 
here for definiteness, although a more general mode ex- 
pansion can be used. We introduce Ai = SM^ fermionic 
or bosonic annihilation operators a — (oks) in momen- 
tum space, labelled by momentum (k = A/cj) and spin 
(s). Here we assume periodic boundaries in a finite vol- 
ume V — , and a lattice of cells, with momentum 
spacing of Afc 271 /L in each coordinate. Localized an- 
nihilation and creation operators a„ on a spatial lattice 



of position indices Tn, with cell volume AV 
defined using a discrete Fourier transform: 



V/M^ 



1 



M3/2 



E 



Oks exp 



27rik . n 
AkM 



(2.4) 



The combined index n is a spin-space 4- vector, n = (n, s). 
In the case of bosonic (fermionic) fields, the commutators 
(anticommutators) are defined as: 



4] 



(2.5) 



The corresponding local number operator is N„i — 
am- The continuum Hamiltonian is regained in the 
limit of a large number of lattice sites of the resulting 
Hubbard type model Hamiltonian: 



H(a\a)^ lim hy^ 



(2.6) 



0J„ 



2m. 



exp 



27rik . (n - n') 



AkM 



5ss' (2.7) 



and the interaction matrix Xmn is (approximately): 

While this is introduced here as an approximation to a 
continuum system, it is also experimentally feasible to 
use an optical lattice (TTI . [l^ to engineer the hamilto- 
nian directly, so that each spatial index coincides with 
a local trapping potential well. In this way one can ob- 
tain a physical system directly corresponding to the fa- 
mous Hubbard model of condensed matter [l3|, except 
that it does not involve the numerous approximations 
that would be required in a typical condensed matter 
setting. In the following calculations, we will assume 
for simplicity that the interaction is local in space, with 
Xnn' — Xss'S^n,; although this is not essential. 



A. Exponential complexity 

The central issue that makes quantum dynamical cal- 
culations difficult in many body quantum physics is the 
issue of exponential complexity |14| . For example, if 
we consider N bosons distributed among Ai modes, the 
number of distinct orthogonal quantum states is obtained 
by combinatorics: how many ways can we divide the 
particles amongst the modes? The number of quantum 
states is then simply: 



(M+N 

{M - ly.Ni 



(2. 



To give relevant numbers, suppose we consider num- 
bers that are approximately typical of many ultra-cold 
atom experiments, with N = Ai = 500,000. One finds 
that: 



i2M 



10 



300,000 



(2.9) 



With fermions, we have fewer states, since each mode can 
have an occupation number of or 1, meaning that: 



10 



150,000 



(2.10) 



In either case, there are more linear equations to solve 
than atoms in the universe. By comparison, the number 
of classical equations would be 



Nc^2M^ 10^ 



(2.11) 



While the classical problem is difficult, it is soluble on 
many current digital computers. The quantum problem, 
on the other hand, is nearly impossible to treat. In par- 
ticular, one can't diagonalize the Hamiltonian, which is 



now a 10 



300,000 



10 



300,000 



matrix in the bosonic case. 
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III. EXACT DYNAMICS 

If there are small numbers of modes, one can indeed 
obtain exact eigenvectors. In this case it is possible to 
diagonalize the Hamiltonian, and obtain the eigenvectors 
and eigenvalues for a small number of particles, typically 
in the range 10 — 100. As an example, we consider the 
generation of entanglement through four-mode nonlinear 
dynamics in two- well trap holding a two-species BEC sys- 
tems, in which the nonlinearity enters through S-wave 
scattering interactions ^5|]. The basic interaction that 
generates entanglement in the first place is the nonlinear 
S-wave scattering interaction, which we consider to be an 
interaction between two spin-states in ^'^ Rb. The two spin 
states labeled i — 1,2 are |1) = |F = 1, nip = +1), |2) = 
\F = 2, mp = —1), and there are two spatial modes cor- 
responding to optical trapping of modes labeled a, b for 
clarity. Thus, S = 2, M = 2 and M ^ SM ^ A. With 
such a small number of modes there is no exponential 
complexity issue, and the Hamiltonian can be exactly di- 
agonalized using numerical techniques. 

The Hamiltonian for the coupled system is: 



H/h- 



+ {a^ o b,} 



(3.1) 

Here uj is the inter-well tunneling rate between the two 
wells with localized modes , bi while Xij is the intra- 
well interaction matrix between the different spin com- 
ponents. 

We can solve this using either Schroedinger or Heisen- 
berg equations of motion. To illustrate this, suppose that 
a; = 0, and we have just one well. In the Heisenberg case, 
one obtains: 



= -iY^XijNja, 



(3.2) 



Since the number of particles is conserved in each 
mode, this has the solution: 



ai (t) = exp 



a^ (0) 



(3.3) 



More generally, it is convenient to use a matrix ex- 
pansion of the Hamiltonian in a number-state basis. For 
dynamics, we explicitly assume that ai, bi and 02, 62 
are initially in coherent states. This models the relative 
coherence between the wells obtained with a low inter- 
well potential barrier, together with an overall Poissonian 
number fluctuation as typically found in an experimen- 
tal BEC. We note that the coherent state also includes 
an overall phase coherence, which has no effect on our 
results. For simplicity, we suppose that the initial state 
is prepared in an overall four-mode coherent state using 
a Rabi rotation: |-0 >= \a >ai |a >bi |a >a2 l^t >b2- 



Next, we assume that the inter- well potential is in- 
creased so that each well evolves independently. Finally, 
we decrease the inter-well potential for a short time, so 
that it acts as a controllable, non-adiabatic beam-splitter 
^16] . to allow interference between the wells, followed by 
independent spin measurements in each well. 



A. Squeezed and entangled states produced by 
double-well BEC 

We now use the techniques given above to investi- 
gate a particular dynamical strategy for generating EPR 
entanglement. The technique treated here is generally 
along the lines investigated experimentally in fibre-optics, 
by comparison with squeezing and entanglement experi- 
ments in optical fibers (l7H19j . An important difference 
is that the fiber experiments use time-delayed pulses to 
eliminate interactions between the components. This is 
not readily feasible in BEC experiments, although Fesh- 
bach resonances can achieve this to some extent. 

Let ai,a2 be operators for two internal states in the A 
well and bi , &2 operators for two internal states at the 
B well. Na — a\a2 + oioi and Nb = b\b2 + b\bi are 
the atom number operators of these modes in each well. 
We define Schwinger spin operators at each site for the 
measurement of the EPR paradox and entanglement. We 
define general, phase-rotated spin components according 
to: 



Jy — 



01026 



oJoie*(«^-«i) 



01026 



J 



o|o2 



oioi ) /2 



(3.4) 



il IS 



at A and similar definition at B, while A9 = 62 
the phase shift between mode 1 and mode 2. 

We the select the phase shift to make sure (Jy) 7^ 0. 
The Schwinger spin operators orthogonal to Jy are given 
as J{9) = cos{9)Jz + sin{9)Jx, all of which have the 
property {J {9}) — 0. We define A9 = tt/2 — a where a 
is time dependent, such that (02O1) — |(o20i)|6*". This 
plane contains an infinite family of maximally conjugate 
Schwinger spin operators, generally given by J{9) and 
J {9 -f 7r/2) which obey the uncertainty relation 



AV(0)AV(0 + ^/2) > \{Jy)\yA 



Thus a state which obeys 



A^J{9) < \{Jy)\/2 < A^J{9 + 7t/2) 



(3.5) 



(3.6) 



is a squeezed state, as shown in the Fig. [21 Here, we 
have optimized the phase choice 9 to get the best squeez- 
ing of the Schwinger spin operators by the criterion that 
dA'^J{9)/d9 — 0, hence obtaining 



tg(20) = 2(J„ J,)/(AV,- AV,). 



(3.7) 
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Figure 2. Squeezing of Schwinger spin operators 
10/ogio (A^ Jo/no) (solid), lOlogio (A^ Jg+^/a/n-o) (dashed), 
and shot noise level no — |(Jj;)|/2 (dash-dotted) via BEC 
with number of Rb atoms. Here the parameters correspond 
to Rb atoms at magnetic field B = 9.131G, with scatter- 
ing lengths ai — 100. 4ao, a22 ~ 95.5ao, and ai = 80.8ao. 
ao = 53pm. The coupling constant Xij oc 2cjxaji- Here 
Na = 200, T = xiiNAt. 



Ill this strategy, spin-squeezing at each site can be 
readily obtained, by unitary evolution from an initial 
coherent state, under the local Haniiltonian. Then we 
can obtain sum and difference spins between two sites: 

a2 jAB _ a2( jA_ jB\ ^ ^ a2 jAB _ a2 ( jA , 

J^^y2)j prior to using the beam-splitter - which is 
achieved by a modulation of the inter-well potential, as 
shown in Fig. ^a). 

After using the beam-splitter, entanglement can be de- 
tected via spin measurements using the spin version of 
the Heisenberg-product entanglement criterion [l^ 
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Figure 3. (a) Squeezing of Schwinger spin operators 
SdB- S+ = 10logio[A^iJ^ -Ji')/no] (solid), S- = 
10/ogio [A^(J^^^/2 + J'^+7r/2)/'^o] (dashed), and no = 
{\{Jy) \ + \ {Jy)\)/'2 is shot noise (dash-dotted), (b) Entan- 
glement (Eproduct) based on the criterion (|3.8p by the solid 
curve and Esum in sum criterion p.9p by the dashed curve. 
(This figure has been published in Ref. 15l |) 



E, 



product 



2 /A2 T^B . A 2 jAB 



\{J^)\ + \U^)\ 

or the sum criterion 1201 



< 1 



(3.8) 



A 2 jAB 



\{Jy^)\ + \{Jy'')\ 



< 1, 



(3.9) 



as shown in Fig. [21[b). 

While diagonalizatioii is possible for even larger parti- 
cle numbers, the two-mode approximation becomes less 
and less reliable. For larger numbers of particles the in- 
teraction energy becomes as large cLS cLS the harmonic 
oscillator mode spacing. This means that few mode 
approximations become inapplicable, and the problem 
quickly develops exponentially large numbers of many- 
body states. Methods for treating this more general case 
are given in the next section. 



IV. CLASSICAL PHASE-SPACE 

Early techniques for calculating the dynamics and 
thermal equilibrium states of large quantum systems used 
techniques based on mappings to a classical phase-space 
[2l|. These were used not just in statistical many-body 
theory, but also in applications involving coherence the- 
ory and lasers. The most widely used methods of this 
type are the Wiener representation[22] and the Glauber- 
Sudarshan [2^ P-represeiitations. These differ in 
some technical details. In essence, both can be used for 
mapping second-quantized bosonic fields to a classical 
phase-space. However, the Wigiier representation cor- 
responds to symmetrically-ordered operator mappings, 
while the P-representation corresponds to normal order- 
ing. It is also possible to use an anti-normal ordered 
mapping, called the Husimi Q- function [2^, but this is 
less commonly used for dynamical calculations. When 
used for quantum fields, in a truncation approximation 
described below, these types of phase-space method are 
sometimes called c-field techniques. Phase-space meth- 
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ods for quantum fields were used to p redict quantum 
squeezing in solitons in fiber optics [26|-|30|, with an ex- 
cellent agreement with subsequent experimental tests (sil- 
[33 |. They were later applied to ultra-cold atomic BEC 
dynamics [34|], and have been widely used, especially at 
finite temperatures [s^ [s^ . 



A. Glauber-Sudarshan P-representation 

This approach uses an overcomplete set of coherent 
states |37l |38|. parameterized by a complex vector a: 



I a) — exp [a^ ■ a — a* ■ a/2\ |0) 



(4.1) 



and one then obtains an expansion for the density matrix 
in the form: 



P 



J P(a)Ai (a) cP^'a (4.2) 



where Ai {a) is a diagonal coherent state projection op- 
erator, which is the basis for the expansion, defined as: 



Ai (a) 



-)"> 



(4.3) 



Here we note that the same expansion can be used for 
either fermions of bosons. In the case of fermions, a is a 
Grassmann variable [9], and one must use the bracketed 
minus sign in Eq (j4.3p . This representation generates 
normal-ordered operator products, in the sense that mo- 
ments of P{a) correspond directly to expectation values 
of normally ordered operator products. 

The advantage of this approach is that it maps quan- 
tum states into M complex coordinates, a — p-f ix, and 
hence only has classical complexity. Another advantage 
is that the use of normally-ordered products means that 
there is no UV vacuum divergence in the expectation 
values. However, for many quantum states, including all 
entangled states, the distribution is not positive, and in- 
deed is highly singular. 

In the bosonic case, the operator basis can be written 
in an alternative form 1391, as: 



A.(A) = 



1 + s 



M 



: exp [-2Sa'<5a/ (1 + s)] : , (4.4) 



where 6a = a — a , (5a^ = a* are relative displace- 
ments, and s indicates the operator ordering (40| . Here 
s = 1 is used for normal ordering, as in the case of the 
P-representation, and other orderings are treated in the 
next subsection. This allows us to recognize that the ba- 
sis is just a Gaussian function of the mode operators: an 
exponential of a quadratic function of annihilation and 
creation operators. Just as with similar Gaussian bases 
for ordinary complex functions, such an operator basis 
has more than one possible form, obtained by changing 
the variance. 



B. Wigner-Moyal phase-space 

An even older phase-space method was developed by 
Wigner [2^, who treated thermal equilibrium problems, 
and Moyal [41|, who extended this to a full dynamical 
equivalence with quantum mechanics. This can also be 
written as an expansion over a Gaussian operator ba- 
sis with symmetric ordering mappings, so that = 0. 
This reduces the basis variance, and therefore increases 
the variance of the distribution function. Formally, the 
expansion is written as: 



•2M, 



(4.5) 



This type of distribution generates symmetrically or- 
dered operator products. It maps quantum states into 
complex or real coordinates, and has the advan- 
tage that this mapping gives dynamical equations that 
are the most similar to classical behavior. Historically, 
Moyal first showed how to map quantum operators into 
differential equations, and had a famous correspondence 
with Dirac, who objected to the fact that the distribution 
had no probabilistic interpretation. 

As result, there is a problem for computational imple- 
mentation. One would like to sample the Wigner distri- 
bution probabilistically, but this is not always possible. 
Since the mapping is nonpositive, there is no generally 
efficient and accurate sampling procedure. The repre- 
sentation is also typically UV divergent in three dimen- 
sions, due to vacuum field fluctuations of symmetrically 
ordered moments. 

For the case of an initial coherent state. 



|*o) = |ao> 



(4.6) 



the initial Wigner distribution is Gaussian and posi- 
tive, so that the quantum noise can be readily sampled 
stochastically, with: 



(4.7) 



where 5olq is a Gaussian random complex number, such 
that the only nonvanishing correlations are: 



1 



(4.8) 



For more general states - even as simple as a num- 
ber state - the Wigner distribution exists but has no 
stochastic equivalent. The quantum dynamical manifes- 
tation of this problem is that one obtains a third-order 
Fokker-Planck equation for the Wigner time evolution 
when there are nonlinear Hamiltonian terms. Such an 
equation has no stochastic equivalent, unless truncated 
to give a second order differential equation. In this ap- 
proximation, the resulting equation has a semi-classical 
form, with: 



.dam 
I — — 
dt 



E 



X-n 



(4.9) 
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While apparently classical, this equation includes quan- 
tum noise in the initial conditions, and can simulate non- 
classical entangled states, in an approximation which is 
valid in the limit of large mode occupations 

It was the nonpositivity of the Wigner representation 
that led Feynman to make his famous conjecture[l3| that 
it was not possible to use a classical or digital computer to 
make a probabilistic representation of a quantum system. 
Nevertheless, while not exact, the truncated Wigner ap- 
proach is relatively simple, and has useful applications as 
a practical approximation when mode occupation num- 
bers are large. 

We note that the Husimi Q-function [2^, which corre- 
sponds to antinormal ordering with O — >■ —1, is always 
positive. Yet, paradoxically, this has no direct stochastic 
equivalent either, since the corresponding Fokker-Planck 
equation is not positive-definite. As we show in the 
next section, despite Feynman's remark, there are routes 
to positive representations of quantum systems that do 
have stochastic equivalents, but they involve enlarged, 
non-classical phase-space mappings. 



C. Large-scale two-component atom interferometry 

As an illustration of the Wigner phase-space tech- 
niques, consider interferometry of a two-component *^Rb 
BEC in a harmonic trap, which is performed by many 
experimental groups worldwide. Atom numbers in these 
experiments range from 10^ to 10^, which makes it im- 
possible to simulate the behavior of the system exactly 
in a few mode approximation: these larger traps are in- 
herently multi-mode. A common approach involves the 
propagation of semi-classical Gross-Pitaevskii equations 
[21| . Although these equations provide a good descrip- 
tion of the condensate evolution, they do not account for 
quantum effects in the cloud, and cannot predict vari- 
ances of the observables. A truncated Wigner phase- 
space ap pro ach can be used to obtain more accurate pre- 
dictions |42|-[4i|. 

The effective Hamiltonian in three dimensions is 



can play a significant part in the evolution, can be in- 
cluded into the model by adding a loss term to the master 
equation (45| : 

^i = --^[H.p]+Y.-^'U d^'^^^'lP]^ (4.12) 

where p is the number of the interacting particles in the 
loss process, the vector 1 specifies the number of particles 
with each spin participating in the interaction, and C is 
the operator functional: 

£1^^ [p] = 20fp0'f^^^0'f'^^0'^^'^p~p0^^^^0'^^'\ (4.13) 

The reservoir coupling operators Oj^'' are the distinct p- 
fold products of local field annihilation operators, 

0[^^ = 0[^) (*) = ^>i, (x)vl/,, (x) . . . (x), (4.14) 

describing local collisional losses. 

After a transformation to the Wigner representation 
using functional Wigner correspondences [43] and trun- 
cation of high-order terms, the resulting equations take a 
form similar to that of coupled GPEs with the addition 
of stochastic terms |36i j: 

u 

-r,(x)+^/3g)C.'^)(x,t), (4.15) 

where (f)s (x) is a Wigner c- field corresponding to the op- 
erator field vE's(x). Additionally, Tg is the nonlinear loss 
term and is the damping noise coefficient which are 

both functions of the Wigner fields, while C'f\yi,t) IS a 
corresponding complex, stochastic delta-correlated Gaus- 
sian noise such that: 



ss' 

(4.10) 

where is an annihilation operator for spin s, with the 
space position omitted for brevity, and Xss' is the spin- 
dependent contact interaction strength. The operator 
Kss' is the single-particle Hamiltonian: 



Kss' = \ -—^^+UJs+Vs{-^)/hj 5ss'+^ss'(t), (4.11) 

where m is the atomic mass, Vs is the external trapping 
potential for spin s, ui^ is the internal energy of spin s, 
VLss' is the time-dependent coupling term. Losses, which 



(C*^) (x, i)4^'^*(x', t')) = 5,^5,,,5^ (X - x') 5 {t t') . 

(4.16) 

These stochastic equations can be solved numerically 
using conventional methods on a discrete lattice. The re- 
sulting equations then have a similar form to the general 
lattice Wigner equations, Eq (|4.9p . apart from additional 
loss and noise terms: 

dam . sr^ r , I 1 2 " 

n 

-r™ + EA?cS(i), (4.17) 

Correlations of any order can be extracted, and specif- 
ically, one can obtain the value of squeezing parameter 
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Figure 4. Wigner simulations of squeezing in the vicinity of 
9.1 G Feshbach resonance in *^Rb plotted with a logarithmic 
scale SdB = Wlogio{£,^). Inter-component scattering lengths: 
ai = 80.0 ao (blue solid line), ai — 85.0 ao (red dashed line), 
ai — 90.0 ao (green dash-dotted line) and ai = 95.0 ao (black 
dotted line). 



This serves as an indicator of entanglement in the 
condensate l46l.l47l: 



NASI, 



(4.18) 



where N is the number of atoms, S is the total spin, and 
AS'^j„ is the minimal variance of spin over all possible 
directions. The squeezing parameter is a fourth-order 
field correlation, with values < 1 indicating entangled, 
or spin squeezed states. By simulating the evolution of 
the condensate under different conditions one can find 
the optimal regime for producing maximum squeezing. 

For example, the inter-component scattering length 
of optically trapped two-component ^^Rb condensate of 
states \F = 1, mp — +1) and \F = 2, iriF = —1) near a 
Feshbach resonance at 9.1 G can be changed by varying 
the strength of magnetic field. When one moves closer to 
the resonance, the inter-component scattering length de- 
creases, providing better squeezing, but inter-component 
losses increase correspondingly [48], eliminating the co- 
herence faster. Fig ^ shows the evolution of squeezing 
parameter in time for four different values of ai, where 
best results are achieved for oi = 90 cq. In other words, 
one has to pick an optimal but non-zero detuning from 
the Feshbach resonance in order to get maximum squeez- 
ing. 



V. NON-CLASSICAL PHASE-SPACE 

We generically regard any expansion of the density ma- 
trix in a complete, non-orthogonal basis set with contin- 
uous variables as a phase-space method. Such methods 
are not restricted to classical mappings, however. In 
the following sections, we review progress in develop- 
ing non-classical phase-space mappings, typically using 
higher than classical dimensionality. This approach is 
essentially a middle ground between the low complex- 



ity of classical phase-space, and the exponentially large 
complexity of the full many-body Hilbert space. 



A. Positive-P function methods 

In this approach, one extends the mapping of a 
bosonic field theory onto classical phase-space used in the 
Glauber-Sudarshan P-function, to a larger phase-space of 
double the classical dimension. This is best thought of 
as a minimal prescription for including coherent state su- 
perpositions and entanglement into the basis set. Thus, 
one defines: 



p = j P+{a.,p)k+{a,l3)(f^'^ad?^''p (5.1) 



where the basis set is now: 



A+ (a,^) = 



(5.2) 



This enlarged phase-space allows positive probabilities 
for any quantum state, since it is possible to prove an 
existence theorem that any physical density matrix has 
a positive distribution in this form. While this itself is 
no different to the properties of the Husimi Q-function 
- another positive representation - there are additional 
advantages, as we explain below. In comparison to the 
usual diagonal Glauber-Sudarshan case, we note that the 
+P representation has these differences: 

• It now maps quantum states into 4M real coordi- 
nates: a, /3 = p + ix, p' + ix' 

• The corresponding phase-space has double the di- 
mensionality of a classical phase-space 

• The advantage is that one can represent superposi- 
tions including entangled states without singulari- 
ties. 



B. +P Existence Theorem 

The most significant property of the +P method is 
the existence theorem: a positive P-function always ex- 
ists, for any density matrix. While the proof is too 
lengthy to be given here, we quote the result, which 
has a simple constructive form [49]. For any hermitian, 
positive-definite density matrix p. there is a correspond- 
ing positive-P distribution, of 'canonical' form: 



P+(a,/3) 



1 



47r2 



M 



-ia-rr/4 



(3* 



(3* 



(5-3) 

The advantage here is not just that the distribution 
is nonsingular, but more importantly that probabilistic 
sampling is possible. This is a crucial issue when dealing 
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with the complexity of many-body systems. Generally, 
random, probabilistic sampling is the only practical ap- 
proach for subduing the unruly nature of an exponen- 
tially complex Hilbert space. 

This approach nevertheless is not without its own prob- 
lems. The use of a non-orthogonal basis means that the 
distribution is non-unique. We can choose a compact 
form - like the 'canonical' positive form given above - 
as an initial condition. However, this form is generally 
not conserved by the equations of motion, since it is not 
unique. If the equations of motion generate distributions 
that are less compact as time evolves, the this allows 
sampling errors to grow with time. 

An important application of the +P distribution is the 
calculation of measurable operator moments. In order to 
calculate an operator expectation value, there is a corre- 
spondence between the moments of the +P distribution, 
and the normally ordered operator products. These come 
directly from the fact that coherent state are eigenstates 
of the annihilation operator, and that Tr [A+(q:,/3)] = 1, 
which means that any normally ordered operator prod- 
uct is simply a stochastic average of the phase-space vari- 
ables: 

{al,--- an) = Tr [• • • a„A+(a, I3)al^, . . .] (5.4) 



C. +P time-evolution 

The route to obtaining time-evolution equations is to 
map operator equations into differential equations for the 
P-function. Differentiating the +P projection operator 
gives the following four identities: 




(5.5) 

Since the projector is an analytic function of both a™ 
and /S™, we can obtain alternate identities by replac- 
ing d/da by either d/dax or d/iday. This equivalence 
allows a positive-definite diffusion to be obtained, with 
stochastic evolution. The result of this procedure is that 
our exponential complex quantum problem is now trans- 
formed into a stochastic equation. Thus, for the case of a 
single-component Bose gas with S-wave interactions, one 
obtains the following equations in the simplest case: 



da"^' 



dt 

.dPni 
I — -- 

dt 



X"m^m + ^m^ (0 am(5.6) 



Here (,m (t) are 2M independent real, random Gaus- 
sian noise terms, with correlations given on a discrete 
spatial lattice with cell volume AV^, by: 



1 



Smm'Sii'S {t - t') . 



(5.7) 



Similar techniques can also be used for the case of 
fermions, using the Gaussian representation method, 
with some modifications which arc not treated in detail 
here. In both cases, the essential trade-off is that, as with 
any sampling technique, many parallel trajectories are 
needed to control growing sampling errors. This can be 
modified by changing the choice of basis, and the stochas- 
tic mapping which is not unique. This approach leads to 
the idea of a 'stochastic gauge' [50], which multiplies the 
basis operator A+ by a random weight i7, and improves 
convergence properties by reducing the sampling error. 

We emphasize that while any stochastic method only 
converges for a large number of samples, the sampling 
error is generally a well-controlled numerical error. Like 
the momentum cut-off, the sample-size can be changed 
and the error monitored using well-defined numerical pro- 
cedures. 



D. Time-reversal tests 

We now consider two examples of the application of 
the +P distribution to quantitative simulations of time- 
evolution under our Hubbard-type Hamiltonian. 

By choosing a modified stochastic gauge, it is possible 
to simulate a very large number of bosons, and verify 
the simulation by carrying out a time-reversal test. This 
takes advantage of the fact that unitary evolution is time- 
reversible, so that changing the Hamiltonian sign will 
cause a reverse evolution to occur, which should recreate 
the initial state. Such tests have been carried out with 
up to lO^'^ interacting bosons [Slj. In Fig ([5]), we show a 
time-reversal test carried out for the single mode anhar- 
monic oscillator with an initial coherent state of a = 10, 
and a corresponding mean initial population oi N = 100 
bosons. The quantity graphed is the mean quadrature 
variable, defined as AT = (a -I- a^) /2. The Hamiltonian 
sign was changed at t — 0.5 , resulting in a restoration 
of the initial coherence. 



Xam/3m + \/'^iX (0 



This calculation also demonstrates a subtle feature of 
the +P stochastic method, which is that the phase-space 
distribution is not unique! In fact, a moment's thought 
serves to illustrate that this is a necessary feature of a 
stochastic method that represents unitary evolution in 
quantum mechanics. If the mapping has stochastic be- 
havior during the time-evolution, then it will spread in 
the phase-space as time evolves. Technically, the phase- 
space entropy must increase. 
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Figure 5. Time-reversal test of aniiarmonic oscillator with 
A'^ — 100 initial bosons in a coherent state. The shaded ar- 
eas for r > 0.5 indicate a slowly growing sampling error-bar, 
calculated using the central limit theorem and the sampled 
variance. 




atoms scattered into 
an ~ spherical shell 



Figure 7. Schematic diagram of a BEC collision. A second 
condensate is produced in situ by a Bragg scattering pulse, 
with a relative velocity of 2vq, causing a quantum collision. 
This results in quantum correlated atoms scattering into a 
spherical shell of velocities around the original mean velocity. 




Figure 6. Phase-space distributions, showing non-uniqueness 
after time-reversal at normalized time r = 0.5. The horizontal 
axis is a product of the phase-space coordinate a and the 
random stochastic weight Q. 



Yet reversing time simply results in another type of 
unitary evolution, also with a stochastic equivalent. The 
result is that the phase-space distribution spreads even 
more, increasing the phase-space entropy yet further. 
This is illustrated in Fig which graphs the distri- 
bution underlying the mean quadratures depicted in Fig 
([5]). Clearly, the final distribution after time-reversal is 
totally different from the initial distribution, which is a 
delta-function in phase-space. The final distribution af- 
ter time- reversal is a Gaussian convolution of the original 
delta-function. Despite this, the two distributions are 
physically identical, and have identical observable mo- 
ments due to the non- uniqueness of the basis set. 



E. BEC collision: lO' bosons, 10^ spatial modes 

Next we consider of extremely large complexity: 

a collision of two Bose condensates, each with a very large 
number of particles and modes. 

Recent experiments on ultra-cold Bose-Einstein con- 
densates have been able to generate collisions of quantum 
condensates with large numbers (> 10^) of particlesjj|. 
These experiments typically use metastable ^iJe* con- 
densates so that he particle correlations can be readily 
measured. The initial state is simply a trapped BEC 
in which half the atoms have been accelerated to a high 
relative velocity compared to the other half using op- 
tical Bragg scattering techniques, as shown in Fig ([7]). 
Atoms collide to produce a scattered halo of correlated 
atoms, involving both spontaneous and stimulated emis- 
sion into the scattered modes. These quantitative experi- 
ments provide a rigorous test of the methodology of these 
simulations. 

We consider the collisionjlll of two pure ^^Na BECs, 
with a similar design to a recent experiment at MIT[5s||, 
and more recent experiments in France (53|. A 1.5 x 10^ 
atom condensate is prepared in a cigar-shaped magnetic 
trap with frequencies 20 Hz axially and 80 Hz radially. 
A brief Bragg laser pulse is used to coherently impart a 
velocity of 2vq = 19.64 mm/s to half of the atoms, that 
is much greater than the sound velocity of 3.1 mm/s. At 
this point the trap is turned off so that the wavepackets 
collide freely. 

The coupling constant g depends on the s-wave scat- 
tering length a (2.75nm in the case of ^^Na). We be- 
gin the simulation in the center-of-mass frame at the 
moment the lasers and trap are turned off (t = 0). 



10 




Figure 8. Comparison of truncated Wigner and +P simu- 
lations, for total scattered particle density in velocity space 
at a given final velocity. While the +P results agree with 
experiments, the results show unphysical negative densities 
and 'ghost' scattering events in the approximate truncated 
Wigner approach. In these simulations, space is discretized 
onto a 432 x 50 x 50 lattice with fcx.max = 1-4 x 10''/m and 
fcy.z.max ~ 6.2 X 10^/m, giving 10^ modes in total, with phys- 
ical parameters given in the text. 



The initial wavefunction is modeled as the coherent- 
state mean-field Gross-Pitaevskii (GP) solution of the 
trapped t < condensate, but modulated with a fac- 
tor ^e^'^'^^ + q-^^qk'^ which imparts initial velocities 
Vx = ±vq = -^hko^ jm in the x direction. 

Collisions have been calculated for up to 10^ modes 
and 10^ interacting bosons[55]. This is a clearly expo- 
nential regime, yet the +P technique is definitely appli- 
cable. There are restrictions on the interaction density 
and total time duration possible before the sampling er- 
ror is too large, but useful results are certainly obtain- 
able, as shown in Fig ([8]). We emphasize that sampling 
error, while easily estimated, is not easily reduced at long 
times due to a rapid growth in the distribution variance, 
which eventually makes these simulations impractical. 

By comparison, with the same parameters the trun- 
cated Wigner method clearly fails to produce physically 
sensible results. There is an uncontrolled error leading 
to a completely unphysical depletion of the vacuum at 
large relative velocity, causing apparently negative parti- 
cle densities which of course cannot occur. This depiction 
leads to 'ghost' particles scattering into low- velocity re- 
gions near the original condensate velocity, which do not 
correspond to any real physical events. 

While the +P simulation produces physically sensible 
results up to the time when sampling errors are too large 
to be useful, the truncated Wigner approach is liable to 
generate completely unphysical behavior. This is due to 
the fact that in a three-dimensional quantum field simu- 
lation, there are a diverging number of unoccupied high- 
momentum modes in the limit of high momentum cutoff. 
These contradict the basic high-occupation number ap- 
proximation inherent in the Wigner truncation. 

In summary, the advantages and disadvantages of 
the +P approach are that it treats exponentially large 
Hilbert spaces without mean-field or factorization as- 
sumptions, including either unitary or non-unitary 



damped evolution. No truncation of the equations of mo- 
tion is required , and there is no UV divergence at large 
k-value. However, for unitary evolution, and especially 
for strong interactions, the sampling error grows in time. 
This is intrinsic to the stochastic method, which leads 
to solutions have relatively large tails. From the funda- 
mental existence theorem, there are known solutions that 
are strongly bounded with small sampling errors, but one 
must use different simulation techniques to access these 
solutions. 



VI. GAUSSIAN REPRESENTATION 

The Gaussian phase-representation is a more general 
phase-space representation than either the Wigner or 
positive-P method. In fact, it includes these as spe- 
cial cases, and extends the phase-space idea to include 
fermions |56l . [s^ . Here we consider a general number- 
conserving Gaussian operator basis, in which any den- 
sity matrix p is expanded in terms of Gaussian opera- 
tors, A(A), defined as exponentials of quadratic forms. 
In order to define the Gaussian operators, we con- 
sider a bosonic or fermionic quantum field with an M- 



= a 



and 



dimensional set of mode operators a) = 

In the bosonic case, we can define &a 
8a) = — /3 as operator displacements, where in gen- 
eral OL and are independent complex vectors. In the 
fermionic case we set these displacements to zero. The 
annihilation and creation operators satisfy (anti) com- 
mutation relations, with (-I-) for fermions and (— ) for 
bosons: 



ai a 



The expansion of the density matrix is: 



P(A)A(A)dA , 



(6.1) 



(6.2) 



where P(A) is the probability density over the phase- 
space, A is the complex vector parameter of the Gaus- 
sian representation , dX is the integration measure, and 
A(A) is a Gaussian operator, defined as a normally or- 
dered exponential of a quadratic form of annihilation and 
creation operators: 

A(A) = ^^^^Jf'- '^^P h'^aVH : (6.3) 

Here, /x is a complex M x M matrix, so that the repre- 
sentation phase space is A = [a, /3, /£] . Af = Tr [A,i(A)] 
is a normalizing factor, and : : indicates normal ordering. 
The normalizing factor has two forms, for bosons and 
fermions respectively: 



Mb = det [/£] ^ 
Nf = det [21 - /£] 



(6.4) 
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The matrix fi is related to the stochastic Green's func- 
tion n as: 



21b = /£ ^ I 



-T 



(6.5) 



In either case, the stochastic average of n over the dis- 
tribution P is physically a normally ordered many-body 
Green's function, so that: 



(6.6) 



Similar methods to the positive-P approach can then 
be used for calculating time-evolution. These have 
been used very successfully in treating, for example, the 
ground state of the fermionic Hubbard model [58] - a 
problem of much interest in condensed matter physics. 



A. Linear entropy 



calculus, we obtain that the inner product of two un- 
normalized fermionic Gaussian operators is: 



F (/£,«.) =det [I + (/-/£) (I-!d)] 



(6.9) 



We rewrite these expressions in terms of the normally 
ordered Green's functions or correlations of the basis sets: 



= Tr 



A(n)a|aj 



(6.10) 



Introducing the hole Green's functions, n = [I — n], and 
m = [I — m] , we obtain the following result for the nor- 
malized inner product of fermionic Gaussian operators: 



Tr [A(m)A(n)] = det [ nm + nm ] 



(6.11) 



Using the result of the trace of the inner product of 
fermionic Gaussian operators, Eq. ()6.1ip . we obtain that 
the expression for the linear entropy in a Gaussian phase 
representation Eq. (|6.8p is: 



Entropy is a measure of loss of information and entan- 
glement in a quantum system, and it is important to have 
a method of sampling a phase-space distribution in order 
to estimate the entropy. Here we show that the Gaussian 
phase-space method is well-suited to this type of calcu- 
lation. In particular, the linear or Renyi entropy ^59|], is 
defined as: 

52 = -lnrr(p2) . (6.7) 

This has similar properties to the usual logarithmic en- 
tropy, and measures state purity, since ^2 = for a pure 
state, while 5*2 > for a mixed state. The Renyi entropy 
in a phase-space representation can be written using Eq. 
(|6.7p and the expansion of the density matrix Eq. (|6.2p 
as: 



S2^-\n jj P{X)P{X')Tr (A(A)A(A')) dXdX' . (6.8) 

In order to obtain an expression of the linear entropy 
using the Gaussian phase-representation is necessary to 
evaluate the inn er p roducts of Gaussian operators of form 
Tr (A(A)A(A'))|60|- We emphasize here that this proce- 
dure does not give useful results for the usual classical 
phase-space representations. However, we will show that 
it yields remarkably simple results for both fermionic and 
bosonic Gaussian representations. 

For the fermionic case, we first evaluate the trace of the 
inner product of two un-normalized fermionic operators 
P (Mi IL) — Tr [A„ (/x) A„ (i/)] . In this approach, the 
physical many-body system is treated as a distribution 
over fermionic Green's functions, whose average are the 
observed Green's function or correlation function. 

In order to evaluate the trace we use fermionic coher- 
ent states \a.) in terms of Grassmann variables Q;[6li, as 
well as the trace of an operator and the identity opera- 
tor of fermionic coherent states. After some Grassmann 



S2 = -lnJJ P{m)P{n) det [nm + nm] dmdn. (6.12) 

Just as in the case of fermions, we evaluate the inner 
product of two un-normalized bosonic Gaussian opera- 
tors, B {n^u) — Tr [A„ (/x) A„ (iv)] , and after using co- 
herent state expansions, we obtain that the inner product 
of two un-normalized bosonic Gaussian operators is: 

B (/£,i/) = det [I -(/£-/) {iL-DY^ (6.13) 

Similar to the fermionic case, we can rewrite this ex- 
pression in terms of the stochastic Green's functions, and 
finally we have that: 



S2=\n jj P{m)P{n) det [I + n + m] dmdn . (6.14) 

In summary, using the results of the inner products of 
Gaussian operators, we can obtain an expression for the 
linear entropy. Since this is just an average over two inde- 
pendent probabilities, it is readily calculable using sam- 
pled phase-space representations. This expression can 
also be used to evaluate the entanglement of a quantum 
system with a reservoir or other coupled system. 

VII. VARIATIONAL METHODS 

While the previous phase-space methods have had a 
long history in physics, there is a different approach with 
an equally long history, namely the use of variational 
techniques. This is a very simple concept, which is that 
one should use an evolution equation which minimizes 
the error. 

The idea was first proposed by Frenkel[62, p436] and 
Dirac. who developed early variational methods. In 
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Dirac's original approach, an effective action method was 
proposed for a variational wave- function \ip (t)), of the 
form: 

6r^S J dt{^ (t)| [ihdt - H] Itp (t)) = (7.1) 

This generates a variational Schroedinger equation, 
which gives the exact Schroedinger equation in the case 
that \ip (t)) is a complete set of wave-functions. In the 
usual application of the method, IV' (^)) is chosen as a 
specific functional form described by a small number of 
free parameters. The properties of this method are that 
results depend on the chosen function, and there are a 
small number of equations. However, one can't easily 
determine errors, and of course, the method doesn't con- 
verge if lip (t)) is incomplete. 

A. Multiconfigurational methods 




Figure 9. What does the multiverse look hke? Here we give 
an illustration of a quantum 'universe' as a superposition of 
coherent states, illustrated by the blue circles. Each blue cir- 
cle represents a different multi-mode coherent state a'"' (t), 
which has an intrinsic uncertainty. The whole quantum state 
or 'universe' is a superposition of M coherent states with dif- 
ferent phases and amplitudes. 



More recent applications of the variational 
approach[63] have focused on the concept of choos- 
ing an expansion for the variational wavefunction that 
is complete in some limit. These are typically sums of 
individual many-body wavefunctions known as configu- 
rations, hence the term multi-configurational approach. 
A common approach in the BEC case is to construct the 
variational wavefunction from sums of multiply-occupied 
Fock states of the form [6^: 

\i;it))... = J2Cnit)ar...alr\0), (7.2) 

n 

where the operators a] are time-dependent operators 
such that: 




Compared to the usual variational approach, the strat- 
egy used here is to systematically increase the Hilbert 
space dimension by increasing the number of modes, with 
full convergence expected as M — > oo . An obvious draw- 
back is that, since M is the number of modes, and the 
Hilbert space is essentially a standard Fock space, there 
is clearly a potential problem with this strategy. There is 
an exponential complex Hilbert space dimension as M in- 
creases. The result of this problem is that the technique 
is currently restricted to one dimension with less than 
100 particles, and relatively weak interactions. However, 
given this restriction, relatively long interaction times are 
possible. 

B. Combining phase-space and variational methods 

Given the success of phase-space approach in dealing 
with complexity, an obvious question is: how can we 



unify variational and coherent states phase-space meth- 
ods? Such an approach would have several potential ad- 
vantages. Coherent states provide description of long- 
range coherence, and in principle are a complete basis. 
The usual stochastic method with finite computational 
samples can cause a large sampling error. However, the 
variational approach can be used to minimize the error. 
The goal of such an approach is to combine a high degree 
of complexity with relatively long time-scale evolution. 

In order to describe the simplest resulting method 
of this type, we recall the Wheeler-Everett multiverse 
concept |65|. That is, suppose we let jo;) be a multi-mode 
coherent state, then an Af component superposition of 
coherent states can be considered a 'multiverse' - a su- 
perposition of Af classical worlds. This is parameterized 
by a coherent 'super- vector X = [<y.^^\ . . . a'-''^'), where 
a — {ao, . . . aAi) is a multi-mode coherent amplitude, 
flo = 1 is the unit operator, and ao is introduced here as 
a combined relative phase and weight parameter, includ- 
ing the state normalization factor. The corresponding 
quantum state defined as: 

TV 
n=l 

This is a constrained multiverse, with a fixed number of 
copies, shown schematically in Fig One cannot keep 
track of all quantum universes! It is also, in quantum 
mechanical terms, a superposition of non-orthogonal co- 
herent states. We find that the use of a variational princi- 
ple introduces interactions between coherent amplitudes. 
The advantage is a much lower sampling error compared 
to independent stochastic evolution. 
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C. Phase-space variational equations 

Having introduced the concept, we now wish to derive 
the resuhing equations [66]. Suppose the Hamiltonian is 



H (a^ , a) , with an exact wavefunction 
interval At — t — to, then clearly: 



ijj ) , and a time- 



-iHAt 



(7.5) 



We wish to minimize local propagation error in calculat- 
ing the computational result]^ (t)), so that: 



6£ = 6\\\i;{t))-e 



-iHAt 



\^{to))\\ =0. 



(7.6) 



Expanding to first order and taking the limit of At 
t — to — >■ 0, gives: 



^{Sip\ [dt+iH] IV') = 0. 



(7.7) 



We introduce \Aip) — \tp (t)) — \ip (to)) and linearize in 
the coherent state parameters, so that: 



|A^) = AX, 



d\ 



(7.8) 



The variational principle then leads to a diff^erential equa- 
tion for the coherent parameters X, of form: 



(7.9) 



This means that there are now P = N{M + 1) differ- 
ential equations to solve, with a combined index /i, where 
the main terms are 



• The variational matrix: Vni, 



The H- vector: = |^ if jV' (t)) . 



1^) 



However, a numerical problem must be treated, which is 
that V^i/ is generally not invertible, owing to the exis- 
tence of multiple minima. We can solve this iteratively, 
in terms of the parameter change AXl^"!, where we set 
AX[°1 = initially and iterate until a stable solution 
is reached. This is similar to the Tikhonov variational 
method[6^]. In detail, suppose we have a set of varia- 
tional parameters X (to), at time t = to. We evaluate the 
variational matrix and Hamiltonian vector a midpoint 
to At/2, by iterating so that XIp-^I = X (to) + AXIp-^I, 
and setting the change in X to AX^pI , where: 



-iAfH[P"^V2- vIp-^IAX 



[p-i] 



(7.10) 



The final step is to propagate to to + At by setting: 
X (to -f At) = X (to) + 2AXW , (7.11) 



in order to move to the next step in time. This method 
leads to stable equations, with no inversion problems. 

In the case of a coherent state expansion, we define an 
energy matrix: 



a reduced amplitude: 



in) 



- (m) 

a) 



= ^kO + [1 - Sko] "fe 

and an inner product: 



(m) 



= exp 



a, 



(m) + 



a, 



(m) 



E 

A:>0 



(m)* (n) 



(7.12) 



(7.13) 



(7.14) 



The variational matrix definitions are then: 



(mn) 



kl 



H 



(m) 



[[1 - 6ko] 6ki + (5, 



~ (m)* ~ (n) 



P 



(mn) 



E 



da 



(mn) 



(7.15) 



For example, in the case of a single variational term, 
one finds equations equivalent to the Gross-Pitaevskii 
mean-field equation, with an additional phase evolution 
for ao: 



dtao = -i[H -a^ ■ dH/da''] 
dta = -idH/da^ . (7.16) 

For a linear Hamiltonian: 



H = sJu:a 



(7.17) 



one finds that each coherent term evolves independently 
of each other term, giving: 



dta^^^ = 



,(") 



(7.18) 



In this case each linear 'universe' is decoupled from the 
others; this is an exact result, and no approximations are 
required. 



D. Recurrences in the anharmonic oscillator 

Finally, we consider the case of a quantum anharmonic 
oscillator, which describes local S-wave scattering inter- 
actions at a single lattice site, so that: 



(7.19) 



This is an extremely strong test of coherent state 
expansion methods. From an initial coherent state 
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Figure 10. Anharmonic recurrence, x quadrature. Blue 
dashed line has 8 components, green solid line has 16 com- 
ponents. This is almost indistinguishable from the the exact 
solution, which is a red dashed line. 



the quantum evolution generates Schroedinger cat su- 
perpositions, with complete recurrences known to oc- 
cur analytically [6^1, as well as experimentally [69j . 
Variational results showing almost error-free dynamics 
throughout the Schroedinger-cat regime oi t = tt, and 
a complete recurrence at t — 2tt are shown in the fig- 
ures. For these numerical results, we used A = 10"-* to 
control the matrix inversion, with four iterations of the 
variational equations at each time point. 

Here we define quadrature variables, 

X = {a + a^)/2 (7.20) 

and 

y = (a-a^)/(2i) (7.21) 

The exact result, given an initial coherent state, is 
known to have recurrences in both quadratures, as shown 
in the graphs. 

While it is simple enough to be analytically soluble, 
this Hamiltonian would lead to large errors with stochas- 
tic phase-space techniques over time-scales comparable 
to the recurrence time. As we show in Figs ([TU|) and 
PTjl . the use of a variational method allows us to track 
the full recurrence, on timescales where any previous 
stochastic phase-space method would give very large er- 
rors [sil, [TQ, [zII ■ Variational convergence is rapid, with 
only small errors visible using Af = 8 components. These 



are almost completely eliminated by using M = 16 coher- 
ent components. This illustrative example is very simple, 
and indeed the case treated here can be solved exactly 
using analytic techniques. However, it does illustrate the 
utility of variational methods in reducing sampling error 
in phase-space simulations. 




Figure 11. Anharmonic recurrence, y quadrature. Parameters 
as in the previous figure. The blue dashed line with small 
oscillations is not quite converged, while the solid green line 
shows excellent convergence to the exact result. 

VIII. OUTLOOK AND SUMMARY 

In summary, there are a growing number of experi- 
ments in ultra-cold atomic physics that probe the world 
of quantum dynamics. Simple, exact results are only pos- 
sible with a small number of interacting modes. Larger, 
complex quantum systems require new techniques to han- 
dle the issue of exponential complexity. While phase- 
space representations using non-orthogonal basis sets can 
can treat highly complex problems, there is often a prob- 
lem with truncations (in the Wigner case) or sampling 
errors that grow in time (in the positive-P case). 

We have shown that there are newer techniques that 
have much promise. Gaussian phase-space methods are 
much more general, and can treat new issues like fermions 
and entropy calculations. Finally, we have derived a 
Tikhonov-based variational approach which is shown to 
dramatically reduce sampling errors in a case which is 
known to be challenging to handle with phase-space 
methods. This approach gives exact results for linear 
couplings. It has lower sampling errors than the stochas- 
tic +P method, as well as reduced variational complexity 
issues compared to multi-configurational approaches in a 
Fock space basis. While the techniques given here are 
only preliminary, this hybrid variational and phase-space 
approach appears promising for future developments. 
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